#pragma once

#include <cstdlib>
#include <cmath>
#include <iostream>
#include <ctime>

using namespace std;

#ifndef		PRINT_LEVEL_WIDTH
#define		PRINT_LEVEL_WIDTH		4
#endif

#ifndef		PI
#define		PI						3.1415926535897932384626f
#endif

#ifndef		EPS
#define		EPS						1e-8f
#endif

#ifndef		MAXIMUM_DOUBLE_NUM_V
#define		MAXIMUM_DOUBLE_NUM_V	1.0e+305
#endif

#ifndef		MAXIMUM_DOUBLE_NUM_C
#define		MAXIMUM_DOUBLE_NUM_C	0.9e+305
#endif

#ifndef		MINIMUM_DOUBLE_EPS_V
#define		MINIMUM_DOUBLE_EPS_V	1.0e-305
#endif

#ifndef		MINIMUM_DOUBLE_EPS_C
#define		MINIMUM_DOUBLE_EPS_C	1.1e-305
#endif

#ifndef		RANDOM_RANGE
#define		RANDOM_RANGE			10000
#endif

namespace loirey
{
	class CMyMath
	{
	public:
		static const double MAXIMUM_LOG_RESULT_C;
		static const double MINIMUM_LOG_RESULT_C;

		template<typename T>
		static T TranslateNumber_Round(double num, T NO_USE)
		{
			return (T)num;
		}
		static int TranslateNumber_Round(double num, int NO_USE)
		{
			return (int)(num<-EPS ? num-0.5 : (num>+EPS ? num+0.5 : 0.0));
		}
		template<typename T>
		static int TranslateNumber_Floor(T num, double eps = EPS)
		{
			return (int)(num+eps) + (num+eps>0 ? 0 : -1);
		}
		template<typename T>
		static int TranslateNumber_Ceil(T num, double eps = EPS)
		{
			return (int)(num-eps) + (num-eps>0 ? +1 : 0);
		}

		static inline double my_exp(double exponent)
		{
			if (exponent > MAXIMUM_LOG_RESULT_C) return MAXIMUM_DOUBLE_NUM_V;
			if (exponent < MINIMUM_LOG_RESULT_C) return MINIMUM_DOUBLE_EPS_V;
			return exp(exponent);
		}
		static inline double my_log(double num)
		{
			if (num < +MINIMUM_DOUBLE_EPS_C)
				return -MAXIMUM_DOUBLE_NUM_V;
			else
				return log(num);
		}
		static inline double log_Product(double log_A, double log_B)
		{
			int aa = ((!(log_A < +MAXIMUM_DOUBLE_NUM_C)) ? +1 : 
				((!(log_A > -MAXIMUM_DOUBLE_NUM_C)) ? -1 : 0));
			int bb = ((!(log_B < +MAXIMUM_DOUBLE_NUM_C)) ? +1 : 
				((!(log_B > -MAXIMUM_DOUBLE_NUM_C)) ? -1 : 0));
			int tt = aa + bb;
			if (tt > 0)
				return +MAXIMUM_DOUBLE_NUM_V;
			else if (tt < 0)
				return -MAXIMUM_DOUBLE_NUM_V;
			else if (!aa && !bb)
				return log_A + log_B;
			else
				return 0.0;
			//if (!(log_A < +MAXIMUM_DOUBLE_NUM_C)) return +MAXIMUM_DOUBLE_NUM_V;
			//if (!(log_B < +MAXIMUM_DOUBLE_NUM_C)) return +MAXIMUM_DOUBLE_NUM_V;
			//if (!(log_A > -MAXIMUM_DOUBLE_NUM_C)) return -MAXIMUM_DOUBLE_NUM_V;
			//if (!(log_B > -MAXIMUM_DOUBLE_NUM_C)) return -MAXIMUM_DOUBLE_NUM_V;
			//return log_A + log_B;
		}
		static inline double log_Power(double log_A, double exponent)
		{
			if (exponent>-MINIMUM_DOUBLE_EPS_C && exponent<+MINIMUM_DOUBLE_EPS_C) return 0.0;
			if (log_A > +MAXIMUM_DOUBLE_NUM_C) return MAXIMUM_DOUBLE_NUM_V;
			if (log_A < -MAXIMUM_DOUBLE_NUM_C) return MAXIMUM_DOUBLE_NUM_V;
			return log_A * exponent;
		}
		static inline double log_Sum(double log_A, double log_B)
		{
			if (!(log_A < +MAXIMUM_DOUBLE_NUM_C)) return +MAXIMUM_DOUBLE_NUM_V;
			if (!(log_B < +MAXIMUM_DOUBLE_NUM_C)) return +MAXIMUM_DOUBLE_NUM_V;
			if (!(log_A > -MAXIMUM_DOUBLE_NUM_C)) return log_B;
			if (!(log_B > -MAXIMUM_DOUBLE_NUM_C)) return log_A;
			double t = log_A;
			if (t < log_B) t = log_B;
			if (log_A - t > MAXIMUM_LOG_RESULT_C) return log_A;
			if (log_B - t > MAXIMUM_LOG_RESULT_C) return log_B;
			double ret = exp(log_A - t) + exp(log_B - t);
			ret = t + log(ret);
			return ret;
		}
		static inline double log_Minus(double log_A, double log_B)
		{
			double ret = log_B - log_A;
			if (ret > -MINIMUM_DOUBLE_EPS_C) return -MAXIMUM_DOUBLE_NUM_V;
			if (ret < -MAXIMUM_DOUBLE_NUM_C) return log_A;

			ret = 1.0 - exp(ret);
			ret = log_A + log(ret);

			return ret;
		}

		static double FormatAngle(double Angle, double Base = 0.0)
		{
			Base -= PI;
			double PI2 = PI+PI;
			if (Angle < Base)
				Angle += PI2 * ceil((Base - Angle) / PI2);
			return Angle - PI2 * floor((Angle - Base) / PI2);
		}

		static double FormatAngle_Small(double Angle, double Base = 0.0)
		{
			Base -= PI/2;
			if (Angle < Base)
				Angle += PI * ceil((Base - Angle) / PI);
			return Angle - PI * floor((Angle - Base) / PI);
		}

		inline static int dcmp(double num)
		{
			return (num<-EPS ? -1 : (num>+EPS));
		}
		template <typename T>	inline static T sqr(const T& x)
		{
			return x*x;
		}
		template <typename T>	inline static T DoInterpolation_Bilinear(T aa, T ab, T ba, T bb, T dx, T dy)
		{
			return aa*(1-dx)*(1-dy) + ab*(1-dx)*dy + ba*dx*(1-dy) + bb*dx*dy;
		}
		template <typename T>	static T Compute_L2_Norm(const T* pData, size_t Len, const T& eps)
		{
			size_t i;
			const T* pd;
			T z;
			z = eps * Len;
			pd = pData;
			for ( i = 0; i < Len; ++i )
				z += sqr(*(pd++));
			z = sqrt(z);
			return z;
		}
		template <typename T>	static void L2_Normalize(T* pData, size_t Len, const T& eps)
		{
			T z = 1.0 / Compute_L2_Norm(pData, Len, eps);
			size_t i;
			T* pd = pData;
			for ( i = 0; i < Len; ++i )
				*(pd++) *= z;
		}

		template <typename T>	inline static T logistic(T x)
		{
			return (T)(1.0/(1.0+my_exp(0.0-x)));
		}

		template <typename T>	inline static T gaussian(T x, T y = 0, T sigma = 1)
		{
			return (sigma<1e-3 ? 1.0 : (exp((x*x+y*y)/(-2.0f*sigma*sigma)) / (sqrt(PI+PI)*sigma)));
		}
		template <typename T>	inline static T gaussian1stDerivative(T x, T sigma = 1)
		{
			sigma = 2.0f*sigma*sigma;
			return -(x+x)/(sigma*sqrt(sigma*PI)) * exp(-x*x/sigma);
		}
		template <typename T>	inline static T gaussian2ndDerivative(T x, T sigma = 1)
		{
			x *= x; sigma = 2.0f*sigma*sigma;
			return (T)((x+x-sigma)/(.5*sigma*sigma*sqrt(sigma*PI)) * exp(-x/sigma));
		}
	};

	class CMyRandom
	{
	public:
		static int _Next()
		{
			int t1 = rand()%512;
			int t2 = rand()%2048;
			int t3 = rand()%2048;
			return ((t1<<22)|(t2<<11)|t3);
		}
	public:
		static void InitRandomSeed()
		{
			srand((unsigned)time(0));
		}
		static void InitRandomSeed(unsigned CertainSeed)
		{
			srand(CertainSeed);
		}
		static int NextInt(int Range)
		{
			if (Range < 1) return 0;
			return _Next()%Range;
		}
		static unsigned __int64 NextInt64(unsigned __int64 Range)
		{
			if (Range < 1) return 0;

			unsigned __int64 ret = 0;
			unsigned __int64 x = _Next();
			unsigned __int64 y = _Next();
			unsigned __int64 z = rand()%2;
			unsigned __int64 M = ((unsigned __int64)32768) * (unsigned __int64)65536;

			ret = z*M*M+y*M+x;
			ret %= Range;

			return ret;
		}
		static double NextDouble()
		{
			double t1 = _Next()/2147483648.0;
			double t2 = _Next()/2147483648.0;
			return t1 + t2 / 2147483648.0;
			//int t1 = rand()%RANDOM_RANGE;
			//int t2 = rand()%RANDOM_RANGE;
			//return 1.0*t1/RANDOM_RANGE + 1.0*t2/RANDOM_RANGE/RANDOM_RANGE;
		}
		static int DiscreteSample(const double* Distribution, int Amount, double SumWeight = 1.0)
		{
			return DiscreteSample(NextDouble(), Distribution, Amount, SumWeight);
		}
		static int DiscreteSample(double rand_pos, const double* Distribution, int Amount, double SumWeight = 1.0)
		{
			double t = rand_pos * SumWeight;
			int i, bad_ret;
			bad_ret = -1;
			for ( i = 0; i < Amount; ++i )
			{
				if (bad_ret<0 || Distribution[bad_ret]<Distribution[i])
					bad_ret = i;
				if (t < Distribution[i])
					return i;
				else
					t -= Distribution[i];
			}
			return bad_ret;
		}
		static void Gaussian2D(
			double xMean, double yMean,
			double xVariance, double yVariance,
			double& xGaussian, double& yGaussian)
		{
			double x1, x2, w;
			for ( w = 2.0; w >= 1.0; )
			{
				x1 = 2.0 * NextDouble() - 1.0;
				x2 = 2.0 * NextDouble() - 1.0;
				w = x1 * x1 + x2 * x2;
			}
			w = sqrt(-2.0 * log(w) / w);
			xGaussian = xMean + x1 * w * xVariance;
			yGaussian = yMean + x2 * w * yVariance;
		}
		static double Gaussian1D(double Mean, double Variance)
		{
			double ret, ty;
			Gaussian2D(Mean, 0.0, Variance, 0.0, ret, ty);
			return ret;
		}
	};

	template <typename TypeOfKernelAndOutput>
	class CKernel1D
	{
	public:
		int	center;
		int	size;
		TypeOfKernelAndOutput* data;
	public:
		CKernel1D(const CKernel1D& b)
		{
			center = b.center;
			size = b.size;
			if (b.data)
			{
				data = new TypeOfKernelAndOutput[size];
				memcpy(data,b.data,sizeof(TypeOfKernelAndOutput)*size);
			}
			else
				data = NULL;
		}
		CKernel1D()
		{
			center = size = -1;
			data = NULL;
		}
		virtual size_t GetMemoryUsage() const
		{
			return sizeof(int) * 2 + sizeof(TypeOfKernelAndOutput*) + sizeof(TypeOfKernelAndOutput) * size;
		}
		void InitBySize(int Size, int Center)
		{
			center = Center;
			size = Size;
			if (data)
				delete[] data;
			data = new TypeOfKernelAndOutput[Size];
		}
		void InitByCenter(int halfSize)
		{
			center = halfSize;
			if (data && (size == center+center+1))
				return;
			size = center+center+1;
			if (data)
				delete[] data;
			data = new TypeOfKernelAndOutput[size];
		}
		void InitBySigma(TypeOfKernelAndOutput sigma)
		{
			if (data)
				delete[] data;
			data = NULL;
			InitByCenter((int)ceil(3.0*sigma+EPS));
		}
		~CKernel1D()
		{
			if (data)
				delete[] data;
		}
	public:
		// (dstData) = (srcData) *convolve* (this kernel)
		template <typename TypeOfInput>
		void _convolve(unsigned int dim, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData, unsigned int dstJumpOffset, unsigned int srcJumpOffset, unsigned int mskJumpOffset) const
		{
			TypeOfKernelAndOutput t;
			TypeOfKernelAndOutput* const pDstHead = dstData;

			const TypeOfKernelAndOutput* const pDstTail = dstData + dim * dstJumpOffset;

			const TypeOfKernelAndOutput* const pMskTail = data + size;
			const TypeOfKernelAndOutput* pMsk;

			const TypeOfInput* const pSrcHead = srcData;
			const TypeOfInput* const pSrcTail = srcData + dim * srcJumpOffset;
			const TypeOfInput* pSrc;

			mskJumpOffset *= srcJumpOffset;
			for (dstData=pDstHead; dstData<pDstTail; dstData+=dstJumpOffset)
			{
				pSrc = srcData - center * mskJumpOffset;
				t = (TypeOfKernelAndOutput)0;
				for (pMsk=data; pMsk<pMskTail; ++pMsk)
				{
					if (pSrc<pSrcHead)
						t += ((TypeOfKernelAndOutput)(*pSrcHead)) * *pMsk;
					else if (pSrc>=pSrcTail)
						t += ((TypeOfKernelAndOutput)(*(pSrcTail-1))) * *pMsk;
					else
						t += ((TypeOfKernelAndOutput)(*pSrc)) * (*pMsk);
					pSrc += mskJumpOffset;
				}
				*dstData = t;
				srcData += srcJumpOffset;
			}
		}
		template <typename TypeOfInput>
		void ConvolveX(unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData, unsigned int dstJumpOffset, unsigned int srcJumpOffset) const
		{
			_convolve(dimX*dimY, dstData, srcData, dstJumpOffset, srcJumpOffset, 1);
		}
		template <typename TypeOfInput>
		void ConvolveX(unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData, unsigned int jumpOffset = 1) const
		{
			_convolve(dimX*dimY, dstData, srcData, jumpOffset, jumpOffset, 1);
		}
		template <typename TypeOfInput>
		void ConvolveY(unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData, unsigned int dstJumpOffset, unsigned int srcJumpOffset) const
		{
			_convolve(dimX*dimY, dstData, srcData, dstJumpOffset, srcJumpOffset, dimX);
		}
		template <typename TypeOfInput>
		void ConvolveY(unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData, unsigned int jumpOffset = 1) const
		{
			_convolve(dimX*dimY, dstData, srcData, jumpOffset, jumpOffset, dimX);
		}
	public:
		void OutputToStream(ostream& outStream, unsigned int stepInto = 0) const
		{
			int i;
			if (stepInto) outStream<<setw(stepInto)<<' ';
			outStream.precision(8);
			outStream<<setw(3)<<size<<':';
			for ( i = 0; i < size; ++i )
			{
				outStream<<' '<<setw(15)<<fixed<<data[i];
				if (i<size-1 && i%4==3)
					outStream<<endl<<setw(stepInto+4)<<' ';
			}
			outStream<<endl;
		}
	public:
		static void GenerateSingle(CKernel1D<TypeOfKernelAndOutput>& DstKernel1D, TypeOfKernelAndOutput scale = (TypeOfKernelAndOutput)1)
		{
			DstKernel1D.InitByCenter(0);
			DstKernel1D.data[0] = scale;
		}
		static void GenerateGaussian(CKernel1D<TypeOfKernelAndOutput>& DstKernel1D, TypeOfKernelAndOutput sigma)
		{
			DstKernel1D.InitBySigma(sigma);
			for ( int i = 0; i <= DstKernel1D.center; ++i )
				DstKernel1D.data[DstKernel1D.center+i] = DstKernel1D.data[DstKernel1D.center-i] =
					CMyMath::gaussian<TypeOfKernelAndOutput>((TypeOfKernelAndOutput)i,0.0f,sigma);
		}
		static void GenerateGaussian1stDerivative(CKernel1D<TypeOfKernelAndOutput>& DstKernel1D, TypeOfKernelAndOutput sigma)
		{
			DstKernel1D.InitBySigma(sigma);
			for ( int i = 0; i <= DstKernel1D.center; ++i )
				DstKernel1D.data[DstKernel1D.center-i] = -(DstKernel1D.data[DstKernel1D.center+i] =
					CMyMath::gaussian1stDerivative<TypeOfKernelAndOutput>((TypeOfKernelAndOutput)i,sigma));
		}
		static void GenerateGaussian2ndDerivative(CKernel1D<TypeOfKernelAndOutput>& DstKernel1D, TypeOfKernelAndOutput sigma)
		{
			DstKernel1D.InitBySigma(sigma);
			for ( int i = 0; i <= DstKernel1D.center; ++i )
				DstKernel1D.data[DstKernel1D.center+i] = DstKernel1D.data[DstKernel1D.center-i] =
				CMyMath::gaussian2ndDerivative<TypeOfKernelAndOutput>((TypeOfKernelAndOutput)i,sigma);
		}
		static void GenerateLinearInterpolation(CKernel1D<TypeOfKernelAndOutput>& DstKernel1D, int Diameter)
		{
			if (Diameter < 1) Diameter = 1;
			double tfv;
			if (Diameter & 1)
			{
				tfv = 1.0 / Diameter;
				DstKernel1D.InitBySize(Diameter, Diameter/2);
				for ( int i = 0; i < Diameter; ++i )
					DstKernel1D.data[i] = tfv;
			}
			else
			{
				tfv = 0.5 / Diameter;
				DstKernel1D.InitBySize(Diameter+1, Diameter/2);
				DstKernel1D.data[0] = DstKernel1D.data[Diameter] = tfv;
				tfv += tfv;
				for ( int i = Diameter-1; i > 0; --i )
					DstKernel1D.data[i] = tfv;
			}
		}
	public:
		void SetSingle(TypeOfKernelAndOutput scale = (T)1)
		{
			GenerateSingle(*this,scale);
		}
		void SetGaussian(TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian(*this,sigma);
		}
		void SetGaussian1stDerivative(TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian1stDerivative(*this,sigma);
		}
		void SetGaussian2ndDerivative(TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian2ndDerivative(*this,sigma);
		}
		void SetLinearInterpolation(int Diameter)
		{
			GenerateLinearInterpolation(*this, Diameter);
		}
	};

	template <typename TypeOfKernelAndOutput>
	class CKernel2D
	{
	public:
		static const int TYPE_CONVOLVE = 1;
		static const int TYPE_PLUS = 0;
	public:
		CKernel1D<TypeOfKernelAndOutput> kernelX;
		CKernel1D<TypeOfKernelAndOutput> kernelY;
		int fType;
		// +1 for convolve, other for plus
	public:
		CKernel2D()
		{
		}
		virtual size_t GetMemoryUsage() const
		{
			return sizeof(int) + kernelX.GetMemoryUsage() + kernelY.GetMemoryUsage();
		}
		void InitBySigma(TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			kernelX.InitBySigma(sigmaX);
			kernelY.InitBySigma(sigmaY);
		}
		void InitBySigma(TypeOfKernelAndOutput sigma)
		{
			InitBySigma(sigma,sigma);
		}
		CKernel2D(TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			InitBySigma(sigmaX,sigmaY);
		}
		CKernel2D(TypeOfKernelAndOutput sigma)
		{
			InitBySigma(sigma,sigma);
		}
	public:
		// fType == +1	: (dstData) = (srcData) *convolveX* kX *convolveY* kY
		// otherwise	: (dstData) =
		//						(srcData) *convolveX* kX *convolveY* kY +
		//						(srcData) *convolveX* kY *convolveY* kX
		template <typename TypeOfInput>
		void Convolve(unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData, unsigned int dstJumpOffset, unsigned int srcJumpOffset) const
		{
			int i, ii;
			ii = dimX*dimY;

			if (fType == +1)
			{
				TypeOfKernelAndOutput* tmp = new TypeOfKernelAndOutput[ii];
				memset(tmp,0,sizeof(TypeOfKernelAndOutput)*ii);
				kernelY.ConvolveY(dimX,dimY,tmp,srcData,1,srcJumpOffset);
				kernelX.ConvolveX(dimX,dimY,dstData,tmp,dstJumpOffset,1);
				delete[] tmp;
			}
			else
			{
				TypeOfKernelAndOutput* tmp1 = new TypeOfKernelAndOutput[ii];
				memset(tmp1,0,sizeof(TypeOfKernelAndOutput)*ii);
				kernelX.ConvolveX(dimX,dimY,dstData,srcData,1,srcJumpOffset);
				kernelY.ConvolveY(dimX,dimY,tmp1,dstData,dstJumpOffset,1);

				TypeOfKernelAndOutput* tmp2 = new TypeOfKernelAndOutput[ii];
				memset(tmp2,0,sizeof(TypeOfKernelAndOutput)*ii);
				kernelY.ConvolveX(dimX,dimY,dstData,srcData,1,srcJumpOffset);
				kernelX.ConvolveY(dimX,dimY,tmp2,dstData,dstJumpOffset,1);

				for ( i = 0; i < ii; ++i )
					dstData[i] = tmp1[i] + tmp2[i];

				delete[] tmp1;
				delete[] tmp2;
			}
		}
		template <typename TypeOfInput>
		void Convolve(unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData, unsigned int jumpOffset = 1) const
		{
			Convolve(dimX,dimY,dstData,srcData,jumpOffset,jumpOffset);
		}
		template <typename TypeOfInput>
		void ConvolveFormatXYZV(
			unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData,
			unsigned int dstDimZ, unsigned int dstZ, unsigned int dstV,
			unsigned int srcDimZ, unsigned int srcZ, unsigned int srcV
			) const
		{
			int xy = dimX*dimY;
			int xyzDst = xy*dstDimZ;
			int xyzSrc = xy*srcDimZ;
			Convolve(
				dimX,dimY,
				dstData+dstZ*xy+dstV*xyzDst,
				srcData+srcZ*xy+srcV*xyzSrc,
				1,1
				);
		}
		template <typename TypeOfInput>
		void ConvolveFormatVXYZ(
			unsigned int dimX, unsigned int dimY, TypeOfKernelAndOutput* dstData, const TypeOfInput* srcData,
			unsigned int dstDimV, unsigned int dstV, unsigned int dstZ,
			unsigned int srcDimV, unsigned int srcV, unsigned int srcZ
			) const
		{
			int xy = dimX*dimY;
			int vxyDst = dstDimV*xy;
			int vxySrc = srcDimV*xy;
			Convolve(
				dimX,dimY,
				dstData+dstZ*vxyDst+dstV,
				srcData+srcZ*vxySrc+srcV,
				dstDimV,srcDimV
				);
		}
	public:
		void OutputToStream(ostream& outStream, unsigned int stepInto = 0) const
		{
			if (stepInto) outStream<<setw(stepInto)<<' ';
			outStream<<'{'<<endl;
			kernelX.OutputToStream(outStream,stepInto+PRINT_LEVEL_WIDTH);
			kernelY.OutputToStream(outStream,stepInto+PRINT_LEVEL_WIDTH);
			if (stepInto) outStream<<setw(stepInto)<<' ';
			outStream<<'}'<<endl;
		}
	public:
		static void GenerateSingle(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput scale = (TypeOfKernelAndOutput)1)
		{
			dstKernel2D.kernelX.SetSingle(scale);
			dstKernel2D.kernelY.SetSingle((TypeOfKernelAndOutput)1);
			dstKernel2D.fType = +1;
		}
		static void GenerateGaussian(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian(dstKernel2D.kernelX, sigmaX);
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian(dstKernel2D.kernelY, sigmaY);
				dstKernel2D.fType = +1;
		}
		static void GenerateGaussian(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian(dstKernel2D,sigma,sigma);
		}
		static void GenerateGaussian1stDerivativeX(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian1stDerivative(dstKernel2D.kernelX, sigmaX);
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian(dstKernel2D.kernelY, sigmaY);
			dstKernel2D.fType = +1;
		}
		static void GenerateGaussian1stDerivativeY(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian(dstKernel2D.kernelX, sigmaX);
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian1stDerivative(dstKernel2D.kernelY, sigmaY);
			dstKernel2D.fType = +1;
		}
		static void GenerateGaussian1stDerivativeX(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian1stDerivativeX(dstKernel2D,sigma,sigma);
		}
		static void GenerateGaussian1stDerivativeY(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian1stDerivativeY(dstKernel2D,sigma,sigma);
		}
		static void GenerateLaplacianOfGaussian(CKernel2D<TypeOfKernelAndOutput>& dstKernel2D, TypeOfKernelAndOutput sigma)
		{
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian(dstKernel2D.kernelX,sigma);
			CKernel1D<TypeOfKernelAndOutput>::GenerateGaussian2ndDerivative(dstKernel2D.kernelY,sigma);
			dstKernel2D.fType = 0;
		}
	public:
		void SetSingle(TypeOfKernelAndOutput scale = (TypeOfKernelAndOutput)1)
		{
			GenerateSingle(*this,scale);
		}
		void SetGaussian(TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			GenerateGaussian(*this,sigmaX,sigmaY);
		}
		void SetGaussian(TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian(*this,sigma,sigma);
		}
		void SetGaussian1stDerivativeX(TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			GenerateGaussian1stDerivativeX(*this,sigmaX,sigmaY);
		}
		void SetGaussian1stDerivativeX(TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian1stDerivativeX(*this,sigma,sigma);
		}
		void SetGaussian1stDerivativeY(TypeOfKernelAndOutput sigmaX, TypeOfKernelAndOutput sigmaY)
		{
			GenerateGaussian1stDerivativeY(*this,sigmaX,sigmaY);
		}
		void SetGaussian1stDerivativeY(TypeOfKernelAndOutput sigma)
		{
			GenerateGaussian1stDerivativeY(*this,sigma,sigma);
		}
		void SetLaplacianOfGaussian(TypeOfKernelAndOutput sigma)
		{
			GenerateLaplacianOfGaussian(*this,sigma);
		}
	};

	//template<typename T>
	//class CRectBox_CenterBased
	//{
	//public:
	//	T cx, cy, w, h;

	//public:
	//	template<typename T2>
	//	void myInit(const T2& cx, const T2& cy, const T2& w, const T2& h)
	//	{
	//		this->cx = CMyMath::TranslateNumber_Round((double)cx, (T2)0);
	//		this->cy = CMyMath::TranslateNumber_Round((double)cy, (T2)0);
	//		this->w = CMyMath::TranslateNumber_Round((double)w, (T2)0);
	//		this->h = CMyMath::TranslateNumber_Round((double)h, (T2)0);
	//	}
	//	template<typename T2>
	//	void myInit(const CRectBox_CenterBased<T2>& AnotherInstance)
	//	{
	//		myInit(AnotherInstance.cx, AnotherInstance.cy, AnotherInstance.w, AnotherInstance.h);
	//	}
	//};

	template<typename T>
	class CRectBox
	{
	public:
		T x,y,w,h;

	public:
		template<typename T2>
		void myInit(const T2& x, const T2& y, const T2& w, const T2& h)
		{
			this->x = CMyMath::TranslateNumber_Round((double)x, (T)0);
			this->y = CMyMath::TranslateNumber_Round((double)y, (T)0);
			this->w = CMyMath::TranslateNumber_Round((double)(x+w), (T)0) - this->x;
			this->h = CMyMath::TranslateNumber_Round((double)(y+h), (T)0) - this->y;
		}
		template<typename T2>
		void myInit(const CRectBox<T2>& AnotherBox)
		{
			myInit(AnotherBox.x, AnotherBox.y, AnotherBox.w, AnotherBox.h);
		}
		template<typename T2>
		CRectBox<T>& operator = (const CRectBox<T2>& AnotherBox)
		{
			myInit(AnotherBox);
			return *this;
		}

		CRectBox()
		{
			x = y = w = h = (T)0;
		}

		template<typename T2>
		CRectBox(const CRectBox<T2>& AnotherBox)
		{
			myInit(AnotherBox);
		}

		//template<typename T2>
		//void TranlateTo_XYLowHigh(int& xL, int& yL, int& xH, int& yH) const
		//{
		//	xL = CMyMath::TranslateNumber_Round((double)x, (int)0);
		//	yL = CMyMath::TranslateNumber_Round((double)y, (int)0);
		//	xH = CMyMath::TranslateNumber_Round((double)(x+w), (int)0);
		//	yH = CMyMath::TranslateNumber_Round((double)(y+h), (int)0);
		//}

		static double CalcAreaOfOverlap(const CRectBox<T>& a, const CRectBox<T>& b)
		{
			double t1 = (a.x>b.x ? a.x : b.x);
			double t2 = (a.x+a.w<b.x+b.w ? a.x+a.w : b.x+b.w);
			double t3 = (a.y>b.y ? a.y : b.y);
			double t4 = (a.y+a.h<b.y+b.h ? a.y+a.h : b.y+b.h);
			if (t2<=t1 || t4<=t3)
				return 0.0;
			else
				return (t2-t1)*(t4-t3);
		}
		static double CalcAO_Percent(const CRectBox<T>& a, const CRectBox<T>& b)
		{
			double s = a.w*a.h + b.w*b.h;
			double t = CalcAreaOfOverlap(a,b);

			return 1.0*t/(s-t);
		}
		double CalcArea() const
		{
			return 1.0*w*h;
		}
		bool InputFromStream(istream& inStream)
		{
			return !!(inStream>>x>>y>>w>>h);
		}
		void OutputToStream(ostream& outStream, bool fPrintEndLine = true) const
		{
			outStream<<setw(4)<<x<<'\t';
			outStream<<setw(4)<<y<<'\t';
			outStream<<setw(4)<<w<<'\t';
			outStream<<setw(4)<<h;
			if (fPrintEndLine)
				outStream<<endl;
			else
				outStream<<'\t';
		}
		void MakeBounded(const T& xLow, const T& yLow, const T& wLow, const T& hLow, const T& xHigh, const T& yHigh)
		{
			if (x < xLow) x = xLow;
			if (x > xHigh-wLow) x = xHigh-wLow;
			if (w < wLow) w = wLow;
			if (w > xHigh-x) w = xHigh-x;
			if (y < yLow) y = yLow;
			if (y > yHigh-wLow) y = yHigh-wLow;
			if (h < hLow) h = hLow;
			if (h > yHigh-y) h = yHigh-y;
		}
	};

	class CSimpleStatistics
	{
	public:
		int Count;
		double SumX;
		double SumXX;

	public:
		CSimpleStatistics()
		{
			Clear();
		}
		virtual void Clear()
		{
			Count = 0;
			SumXX = SumX = 0.0;
		}
		virtual void NewInstance(double X)
		{
			++Count;
			SumX += X;
			SumXX += X * X;
		}
		virtual void Analyse(double& Mean, double& Var)
		{
			if (Count < 1)
			{
				Mean = Var = 0.0;
				return;
			}
			Mean = SumX / Count;
			Var = SumXX / Count;
			Var = Var - Mean * Mean;
			Var = (Var<MINIMUM_DOUBLE_EPS_C ? 0.0 : sqrt(Var));
		}
	};

}

